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Abstract 

We consider a recently discovered representation for the general solu- 
tion of the Sturm-Liouville equation as a spectral parameter power series 
(SPPS). The coefficients of the power series are given in terms of a par- 
ticular solution of the Sturm-Liouville equation with the zero spectral pa- 
rameter. We show that, among other possible applications, this provides a 
new and efficient numerical method for solving initial value and boundary 
value problems. Moreover, due to its convenient form the representation 
lends itself to numerical solution of spectral Sturm-Liouville problems, 
effectively by calculation of the roots of a polynomial. We discuss exam- 
ples of the numerical implementation of the SPPS method and show it 
to be equally applicable to a wide class of singular Sturm-Liouville prob- 
lems as well as to problems with spectral parameter dependent boundary 
conditions. 

Introduction 

In the recent work [M] a representation was obtained for solutions of the equa- 
tion 

{pu'Y + qu = Xu (1) 
in terms of a known non-trivial solution of the equation 

{pu'o)' + quo = 0, (2) 

where p, q, u, uq are complex-valued functions of the real variable x satisfying 
certain smoothness conditions, and A is an arbitrary complex constant. The 
general solution of ([T]) has the form of a power series with respect to the spectral 
parameter A. The fact that under certain conditions it is a complex analytic 
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function of A has long been known (see, e.g., [H]). The representation from [Tl] 
gives a simple procedure for constructing its Taylor coefficients. 

Our main purpose here is to show that this representation gives us, among 
other possible applications, a simple and powerful method for numerical solu- 
tion of initial value, boundary value and spectral problems. We will generalize 
slightly by considering the equation 

{pu)' + qu — Xru (3) 

and representing its general solution in terms of a solution of which does 
not depend on the coefficient r. 

In Section[T]we prove the main theorem and make several observations which 
will apply to the subsequent applications. In Section [2] we observe that the rate 
of convergence of the spectral parameter power series (SPPS) corresponding 
to ([3]) is quite easy to estimate. This fact and the simple form of the so- 
lutions make the SPPS method appropriate for solving initial and boundary 
value problems. Perhaps more importantly, the representation of a solution of 
the Sturm-Liouville equation in the form of an SPPS is well-suited for solving 
spectral problems as it reduces problem to calculation of zeros of an analytic 
function defined by its Taylor series. In Section [3] we show that even difficult 
examples like the Coffey-Evans potential can be handled by the SPPS method; 
there is little doubt that the numerical performance of the method can be sig- 
nificantly improved using additional computational techniques discussed in the 
paper. In Section U] we show that the appearance of the spectral parameter in 
the boundary conditions does not provoke additional difficulties in application 
of the SPPS method. As an example of application of the SPPS method to 
singular Sturm-Liouville problems we consider in Section [5] a very interesting 
example of a highly non-self-adjoint operator (with complex and singular coeffi- 
cients) discussed in a number of recent publications [U [31 [71 [TOl HH US [H] . The 
SPPS method makes it possible to obtain eigenvalues and eigenfunctions of the 
problem even when other known algorithms encounter numerical difficulties. 



1 Solution of the Sturm-Liouville equation 

Here we formulate a generalization of the main result of [14 , which allows us 
to obtain a general solution of the Sturm-Liouville equation ([3]) in the form of a 
spectral parameter power series. The proof does not depend on pscudoanalytic 
function theory (which was the tool employed in p^). 

Theorem 1 Assume that on a finite interval [a,b], equation (0) possesses a 
particular solution uq such that the functions u^r and l/(ugp) are continuous 
on [a, 6]. Then the general solution of fji on (a, 6) has the form 

U = CiMi + C2W2 (4) 
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where ci and C2 are arbitrary complex constants, 

oo oo 

= wo^A'=^(2'=) and U2 = uo^A'^X(2^-+i) 

fc=0 fc=0 

with and X^^^ being defined by the recursive relations 



(5) 



(6) 



X("-i)(s)w2(s)r(s)ds, n odd, 



^Hs)^-rT-rT ds, n even, 



( X 



odd, 



X = < 



X("-i)(s)u^(s)r(s)ds, n 



K. xq 



(7) 



(8) 



where xo is an arbitrary point in [a, b] such that p is continuous at xq and 
p{xo) 7^ 0. Further, both series in (0) converge uniformly on [a,b]. 

Proof. First we prove that ui and U2 are indeed solutions of ([3]) whenever the 
application of the operator L — -j^p-^ + 9 to them makes sense. For this, note 
that if Lun = 0, then L can be written in the factorized form L = —-j-pu?.-^ — . 

^ ' U(j ax ^ ^> ax uo 

Then application of to ui gives 
1 



-Lui 

r ruQ ax 



1 ^ 'p„2 4.y-;^fe^(2fe) 



dx'^ 

k=Q 



ruQ dx ^-^ 



k=l 

In a similar way one can check that U2 satisfies ([3|) as well. In order to give 
sense to this chain of equalities it is sufficient to prove the uniform convergence 
of the series involved in ui and U2 as well as of the series obtained by a term-wise 
differentiation. This can be done with the aid of the Weierstrass M-test. Indeed, 

k " 



we have 



X(2k) 



< 



(max I rug [)'' ( 



max 



convergent where 



^72l)r and the series J^Wy. 

k=0 



c = |A| (max j^'Mol) ^max 



pul 



\b-a\' 



(9) 
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The uniform convergence of the series in U2 as weh as of the series of derivatives 
can be shown similarly. 

The last step is to verify that the Wronskian of ui and 1*2 is different from 
zero at least at one point (which necessarily implies the linear independence of 
Ml and U2 on the whole segment [a, b]). It is easy to see that by definition all the 
X(")(xo) and X(")(xo) vanish except for X^°\xo) and X(°)(cco) which equal 1. 
Thus 

ui{xo) = uo{xo), u[{xa) = Uo(xo), (10) 
U2{xo) = 0, u'^ixo) = \ — - (11) 

Uq{Xq)p{xq) 

and the Wronskian of ui and U2 at xq equals l/p{xo) 7^ 0. ■ 

Remark 2 In the case A = 0, the solution ^ becomes u\ = uq and U2 = 

„„ / . , ... „ . „ 

Xo 

second linearly independent solution. 

Remark 3 The result of Theorem [7] is valid for infinite intervals as well, the 
series being uniformly convergent on any finite subinterval. 

Remark 4 One of the functions rwg or 1/{pUq) may not be continuous on [a, b] 
and yet ui or U2 may make sense. For example, in the case of the Bessel equa- 
tion {xu'y — = —Xxu, we can choose uq{x) = x/2. Then 1/{puq) ^ C[0, 1]. 
Nevertheless all integrals in ^ exist and ui coincides with the nonsingular 
Ji{^/\x), while U2 is a singular solution of the Bessel equation. 

Remark 5 In the regular case the existence and construction of the required 
Uq presents no difficulty. Let p and q be real valued, p{x) 7^ for all x G [a, b] 
and let p, p' , r and q be continuous on [a,b]. Then (0j possesses two linearly 
independent regular solutions vi and V2 whose zeros alternate. Thus one may 
choose Uq = Vi + iV2. 

Even when the coefficient functions are not real-valued, nonvanishing solu- 
tions abound. Assume that the equation possesses a pair of linearly independent 
solutions vi and V2. Let Aj = {x G {a,b): Vj{x) ^ 0}, j = 1,2. The function 
tp{c, x) — cvi(x) -\-v2ix) for c G C, X G Ai, has maximal (real) rank 2 (i.e., ij) is 
a submersion) if and only if vi[x) 7^ 0. By the preimage theorem '13, chapter 
1], the preimage '!/;^^(0) is a 1-dimensional submanifold ofC x Ai and therefore 
its projection i?i C C has zero measure. For every c G C\-Bi, the linear combi- 
nation cvi + V2 does not vanish in Ai . The same reasoning gives us a null set 
B2 such that vi + cv2 does not vanish in A2 when c G C\B2. Altogether, there 
exist nonvanishing solutions throughout Ai U A2, which is all of (a, b) because 
vi, V2 are linearly independent solutions of a second order differential equation. 
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Remark 6 The procedure for construction of solutions described in Theorem]^ 
works not only when a solution is available for A = 0, but in fact when a solution 
of the equation 

(pUgY + quo = Xoruo (12) 
is known for some fixed Xq. The solution (0) now takes the form 



-(2fe+l) 

k=0 fc=0 

This can be easily verified by writing (0) as 

{L — Xor) u = (A — Ao) ru. 

The operator on the left-hand side can be factorized exactly as in the proof of 
the theorem, and the same reasoning carries through. 

Remark 7 For calculating the series in it may be convenient to calculate 
Qj- xi"^) directly from X^""^) or X*^"^^-'. For example, when n is even we 

have 

r^TwT f no{tfr{t)X^^^-^\t)dtds 

Jxo uo{syp[s) j^^ 

= [{Fix) - P{t))uo{tfr{t)X^''-^\t)dt 

where P' = 1/(wqp). 

Remark 8 Other representations of the general solution of fjj as a formal 
power series have been long known (see Theorem 1], ^) and used for 

studying qualitative properties of solutions. The complicated manner in which 
the parameter X appears in those representations makes that form of a general 
solution too difficult for quantitative analysis of spectral and boundary value 
problems. In contrast, the solution is a power series with respect to X, 

making it quite attractive for numerical solution of spectral, initial value and 
boundary value problems. 

A special case of Theorem [H with g = 0, A = 1, was known to H. Weyl. 

Corollary 9 ( |23| ) Let 1/p and r be continuous on [a, h\. The general solution 
of the equation 

{pu'Y = ru (13) 

on (a, b) has the form 

U = CiMi + C2W2 (14) 

where c\ and ci are arbitrary constants and u\, ui are defined by \^-\E^ with 
uo = X — \. 

This coroUary enables us to find the particular solution mq discussed in Re- 
mark [5j 
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2 Numerical solution of initial value problems 



Consider the Sturm-Liouville equation ([3]) on [a, b] with any desired initial con- 
ditions. The numerical implementation of the solution via the representation 
([S]) for a general solution is algorithmically simple. One must consider the ac- 
curacy of calculation of the iterated integrals in ([7]) and and the rate of 
convergence of the series ([5]), because in numerical work one must work with 
finitely many terms. 

The main parameters that one can control are the number M of subintcrvals 
in which to divide [a, b] when integrating numerically and the number N of 
powers in the truncated series. The relationship of M to the accuracy of the 
integrals is a standard question and will not be discussed here. In regards to 
N, observe that one can not always expect a good approximation to u over 
all of [a, b] with a series of N terms, no matter how precisely the integrals are 



and Ix'^^'' -'^'1 (see the proof 



calculated. However, using the estimate for 
of Theorem [T] below) it is easy to obtain a rough but useful estimate for the 

N 

where ui^n = uq^X'^X^'^''^ We 
fe=o 



tail of the SPPS. Namely, consider \ui — ui^n] 
have 



\ui - ui^n\ = \uo\ 



k=N+l 



(2k) 



max \ uo\ 



< max |uo| 

N 

cosh ^/c — 



(2fc)! 

k=N+l^ ' 



k=0 



{2k)\ 



where c is defined by ([5]) . In a similar way one can see that the remainder of the 
SPPS corresponding to U2 is estimated by the tail of the power series of sinh -^/c. 
Thus, if a certain value of N is seen to be insufficient for achieving a required 
accuracy, the interval can be subdivided and the initial value problem solved on 
the first subinterval. The initial values of the solution for the second subinterval 

oo 

are calculated easily taking into account that u'-, — —ui H — —y X'^xi^''-^) 

fc = l 

(and analogously for M2). Thus, no numerical differentiation is necessary and 
this process can be continued with little loss in accuracy. 

The required particular solution uq may be calculated using any available 
algorithm; in the examples presented below we have applied the formula of 
Corollary m applying the above subdivision procedure. All of the calculations 
were performed with Mathematica (Wolfram). 



3 Spectral problems 

The fact that spectral Sturm-Liouville problems are related to the problem 
of finding zeros of complex analytic functions of the variable A is quite well 
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known (see, e.g., [H]). For a regular Sturm-Liouville problem the corresponding 
analytic function is even entire. The representation (|4])-(l8]) allows us to obtain 
the Taylor series of that analytic function explicitly. As an example, let us first 
consider a spectral problem for ([3]) with the boundary conditions 

m(0) = and = 0. (15) 

We suppose that the coefficients satisfy the conditions from Remark [S] and that 
uq is constructed as described there, taking xq — 0. From the first boundary 
condition and PTI)) . the constant ci in Q must be zero. Then the spectral 

problem reduces to finding values of A for which U2(l) = uo(l)^^A''X*^^'^+^^ (1) 

k=0 

vanishes. In other words, this spectral problem reduces to the calculation of 

oo 

zeros of the complex analytic function k(A) = ^^OmA™ where 

Now let a and /3 be arbitrary real numbers and consider the more general 
boundary conditions 

u(a) cos a + u'(a) sin a — (16) 

u(6)cos/3 + w'(6)sin/3 = (17) 

together with equation ([3]). Taking the solutions ui and U2 defined by ([5]) and 
using pH)) . pT|) with xq = a, we obtain from (|16p the following equation. 

Sill Ct 

ci (uq (a) cos a + Uq (a) sin a) + C2 — , . , . ~ 0, 

uo{a)p(a) 

which gives C2 = 7C1 when a ^ im, with 7 — ~uo{a)p{a){uo{a) cot a + Mo(a)), 
whereas ci — when a — im. In the latter case the result is similar to the 
example considered above, thus let us suppose a ^ nn. From the definition of 
ui and U2 we have 



^ + A'=l(2fe-i) and u'^^U2 + —yX''X 

Then the boundary condition (fT7|) implies that 

/ 00 00 ^ 

{uo{b) cos 13 + u'oib) sin/3) "^X^x'-^^^b) + 7^A^-x(2fc+i)(6) 



(2fc) 



^fe=0 fc=0 



r, I 00 00 \ 



"o(fo)p(&) 
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Thus the spectral problem ([3]), ([TT]) reduces to the problem of calculating 

oc 

zeros of the analytic function k(A) = QmA™ where 

m=0 

ao = K(6) cos/3 + u^fe) sin/3) (1 + 7^^'^(&)) + 

Uo{b)p[b) 



and 



: (moW cosp + u'oib) sin/3) + 7X(2™+i)(6) 

sin/3 (^^(2™-i)(5)+^;f(2m)(^)^ ^ m = l,2,.... 



uo(&)p(6) 

This reduction of a Sturm-Liouville spectral problem lends itself to a simple 
numerical implementation. To calculate the first n eigenvalues we consider 

N 

the Taylor polynomial — ^^^imA™ with N > n. Thus the numerical 

m— 

approximation of eigenvalues of the Sturm-Liouville problem reduces to the 
calculation of zeros of the polynomial KAr(A). 

There is no need to work with zeros of only one polynomial. It is well known 
that in general the higher roots of a polynomial become less stable with respect 
to small inaccuracies in coefficients. Our spectral parameter power series method 
is well suited to overcome this problem and thus to calculate higher eigenvalues 
with a good accuracy. This is done using Remark [G] Suppose we have already 
calculated the eigenvalue Ao using the procedure described above as a first root 
of the obtained polynomial. Then for the next step we define Uq — ui + iu2 
where ui and U2 are defined by (O with A = Aq. The function Uq is then a 
solution of l|12p. We use it to obtain the eigenvalue Ai of the original problem 
observing that Ai = Ai + Aq where Ai is the first eigenvalue of the equation 
(L — Xor)u = Au with the same boundary conditions as in the original problem. 
This procedure can be continued for calculating higher eigenvalues. Note that 
if Ao = we should begin this shifting procedure starting with Ai. 

Here we discuss some numerical examples. 

Paine Problem. A number of spectral problems which have become standard 
test cases appear in [18l [19] . As a first example we consider 



p{x) = -1, q{x) 



(X + 0.1)2' 

u(0) = 0, u(7r) = 0. 

The eigenvalues in the following table were calculated via SPPS using in- 
tegration on 10,000 subintervals for calculating N = 100 powers of A. These 
eigenvalues were found as roots of a single polynomial (i.e., the shifting of A as 
described in Remark [5] was not applied) . Due to the sensitivity of the larger 
roots of the polynomial to errors in the coefficients, 100-digit arithmetic was 
used. 
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n 


An .19J 


A„ SPPS 





1.5198658211 


1.519865821099 


1 


4.9433098221 


4.943309822144 


2 


10.284662645 


10.28466264509 


3 


17.559957746 


17.55995774633 


4 


26.782863158 


26.78286315899 


5 


37.964425862 


37.96442587941 


6 


51.113357757 


51.11335707578 


7 


66.236447704 


66.23646092491 


8 


83.338962374 


83.33879073183 


9 


102.42498840 


102.4259718823 


10 


123.49770680 


123.512483827 



On the basis of the above values, a new calculation was made by shifting 
with A* = 66, resulting in the following improved approximations for the last 
few eigenvalues. 



n 


A„ ^ 


A„ SPPS 


7 


66.236447704 


66.23644770359 


8 


83.338962374 


83.33896237419 


9 


102.42498840 


102.42498839828 


10 


123.49770680 


123.49770680101 


11 


146.55960608 


146.55960605783 


12 


171.61264485 


171.61265439928 



With A* — 146 and increasing the number of powers to = 150, the 
following further values were obtained. 



n 


A„ m 


A„ SPPS 


11 


146.55960608 


146.55586199495330 


12 


171.61264485 


171.60875781110985 


13 


198.65837500 


198.65416389844202 



When the number of digits for internal calculations was increased to 150, 
SPPS produced the same results. 

Coffey-Evans equation. This test case, defined by 

p{x) ^ ~l, q{x) ^ -2l3cos2x + P"^ sin^ 2x. 

u(-7r/2) = 0, u(7r/2) = 0, 

presents the challenge of distinguishing eigenvalues within the triple clusters 
which form as the parameter f3 increases. We present results for /? = 20, 30, 50. 
In all cases given here the eigenvalues were obtained without shifting A. 
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P = 20. 

M = 10, 000 subintervals, N= 180 powers, 100 digits of precision. 



n 




A„ SPPS 





-0.00000000000000 


0.0000000000000003 


1 


77.91619567714397 


77.9161956771439703 


2 


151.46277834645663 


151.4627783464566396 


3 


151.46322365765863 


151.4632236576586490 


4 


151.46366898835165 


151.4636689883516575 


5 


220.15422983525995 


220.1542298352599497 


6 


283.0948 


283.0948146954014377 


7 


283.2507 


283.2507437431126800 


8 


283.4087 


283.4087354034293064 



/3 = 30 

M ~ 10, 000 subintervals, N = 150 powers, 100 digits of precision. 



n 


A„ ^15, 19] 


A„ SPPS 





0.00000000000000 


0.000000000000000002 


1 


117.946307662070 


117.94630766206876 


2 




231.664928928423790 


3 


231.66492931296 


231.664928928423791 


4 




231.664930082035462 


5 




340.888299091685489 


6 




403.219684016171863 


7 




403.219684016171917 


(3 = 50 

.0, 000 subintervals, N = 150 powers, 100 digits of pre 


n 


A„ ,19^ 


A„ SPPS 





0.00000000000000 


0.000000000000000003 


1 


197.968726516507 


197.96872651650729 


2 




391.807 


3 


391.80819148905 


391.810 


4 




547.1397060 



4 Sturm-Liouville problems with spectral param- 
eter dependent boundary conditions 

In this section we consider Sturm-Liouville problems of the form 

ipu'y + qu = Xru, X £ [a, 6], (18) 

u(a) cos a + u'(a) sin a = 0, Q!e[0, tt), (19) 

/3iu(6) - H^v^ih) = ^(A) {fi\u{h) - (3'^v!{h)) , (20) 

where is a complex- valued function of the variable A and /3i, /32, /3i, *re 
complex numbers. This kind of problem arises in many physical applications 
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(we refer to [T] and references therein) and has been studied in a considerable 
number of publications [U O [H [9l [121 [20] . For some special forms of the function 
V5 such as (f{X) = A or ip{X) = A^ + ciA + C2, results were obtained [H], [10] 
concerning the regularity of the problem (fT51) - ([^ : we will not dwell upon the 
details. Our purpose is to show the applicability of the spectral parameter power 
series (SPPS) method to this type of Sturm-Liouville problems. For simplicity, 
let us suppose that a = and hence the condition ([19]) becomes u{a) = 0. Then 
as was shown in the preceding section, if an eigenfunction exists it necessarily 
coincides with U2 up to a multiplicative constant. 

In this case condition ([^0]) becomes equivalent to the equality 

oo / . , oo 

(uo(6)^i(A) - u',ib)MX)) J2x'x(''^+^\b) - -^^A^X(2'=)(6) = 

(21) 

where (^i2(A) = /3i 2 — (3'i2fW- Calculation of eigenvalues given by pT]) is 
especially simple in the case of being a polynomial of A. Precisely this partic- 
ular situation was considered in all of the abovementioned references concerning 
Sturm-Liouville problems with spectral parameter dependent boundary condi- 
tions. For these problems the calculation of eigenvalues using our method does 
not present any additional difficulty compared to the parameter independent 
situation discussed in the preceding section. 



5 Singular problems 

As was mentioned in Remark[4j one of the solutions of the Sturm-Liouville equa- 
tion can be singular and nevertheless the method presented here is still applica- 
ble. We show one such application to an interesting problem first considered in 
|2] and then in a number of recent pubHcations [H [T] [TO] [TT] [21] [22] . We consider 
on the interval (— tt, tt) the singular non-symmetric differential equation 

. d { . du\ .du ^ 
— i£— smx— — I— — \u (22) 
da; \ dx J dx 

with < e < 2 and periodic conditions at — tt and vr. In spite of the fact that 
([22]) is highly non-self-adjoint (with complex and singular coefficients) all of the 
eigenvalues are real [ST and the spectrum is discrete [lOj. Several algorithms 
have appeared for approximating the eigenvalues to which we compare our re- 
sults. It is known (e.g., |3j) that for each A G C, equation ([22]) possesses a unique 
(up to scalar multiples) solution in L^(— tt, tt) which we denote by ip(x, A). This 
solution may be normalized by the condition (p{Q, A) = 1. Any solution linearly 
independent of Lp will blow up as x — s- 0. Also [3 A is an eigenvalue if and only 
if 

.^(-TT, A) = <^(7r,A). 
It is not difficult to rewrite equation ([22]l in the form ([3]), 
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Choosing uq = 1 as a particular solution corresponding to A = we immediately 
find that ui defined by ([5]) belongs to L^(— tt, tt), while U2 is the singular solution 
blowing up as X ^ 0. Moreover, due to the particular form of the coefficients 
in we have from Remark [7] that 



X(^'=)(a;) = -(tan|) 
fc = l,2,.... 



-l/e 



tan0'^'x2('=-i)(s)(is. 



Both integrals are well behaved at a; = tt. Moreover, the integral (tan |) ds 
can be expressed explicitly in terms of hypergeometric functions (produced for 
example by symbolic manipulations in Mathematica version 6). In the numerical 
calculation of X^^'^'> and subsequently of ui we may approximate the first integral 
above by the sum of corresponding integrals on subintervals, each one calculated 
as follows: 

[\ (tan|)^^^X^(-^)(.)...X^('=-)(:^^) £^ (tani)^^^ ds. 

We show the results of application of the SPPS method in comparison with some 
values calculated previously. For e — 0.5 our values lie between those previously 
published. For e = 0.1 our results display a remarkable agreement with those 
of 

Eigenvalues of ([22)) for e = 0.5 



n 


A„ [22 


A« [f 


A„ SPPS 


1 


1.16714 


1.167342 


1.16723 


2 


2.96821 


2.968852 


2.96844 


3 


5.48168 


5.483680 


5.48268 


4 


8.71272 


8.715534 


8.71354 


5 


12.66119 




12.6618 


6 


17.32643 




17.3275 


7 


22.71033 




22.7110 


8 


28.81106 




28.8122 


9 


35.62928 




35.6311 


10 


43.16666 




43.1677 
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Eigenvalues of ([^ for e = 0.1 



n 


An [41 


An ^^^^ 


An [iUj 


An JrrD 


1 

L 


1 nnQ7 


1 nnQzin 

1 .uuy^u 


i.uuyuo 


i.uuyuo 


2 


2.0733 


2.07305 


2.07334 


2.07334 


3 


3.2297 


3.22894 


3.22978 


3.22978 


4 


4.5012 


4.50088 


4.50134 


4.50134 


5 


5.8992 


5.89968 


5.89993 


5.89993 


6 


7.4298 


7.43154 


7.43194 


7.43194 


7 


9.0951 


9.10034 


9.10097 


9.10097 


8 


10.8945 


10.90881 


10.9092 


10.9092 


9 


12.8252 


12.85742 


12.8578 


12.8578 


10 


14.8820 


14.94727 


14.9478 


14.9478 


15 






27.5331 


27.5331 


20 






43.74 


43.6923 



We also show the results we obtained for e — 0.01, for which we do not know 
any previous result. For this reason we include information on the accuracy 
of the eigenfunctions, estimated as follows. The calculated eigenfunctions u 
were effectively substituted into the differential equation by applying the 
corresponding integral operators which in principle should produce Au, and the 
discrepancy 6i at the right endpoint tt was tabulated. Then the discrepancy 62 
in the right boundary condition was evaluated. 



Eigenvalues of (|^ for e = 0.01 



n 


An SPPS 


Si 


S2 


1 


1.0001 


1.0x10-1-* 


3.8xl0~i« 


2 


2.0008 


1.5x10"^^ 


1.7x10-1'* 


3 


3.00269 


3.0x10-1^ 


7.3x10-1* 


4 


4.00638 


4.9x10-1* 


2.9x10-*^ 


5 


5.01243 


6.8x10-1" 


1.3xl0-*i 


6 


6.02143 


8.8x10-'-* 


6.5x10-'-* 


7 


7.03393 


3.2x10-"^ 


6.7xl0-« 


8 


8.05048 


6.2xl0-« 


3.4x10-'^ 


9 


9.07162 


0.000090 


4.0x10-^ 


10 


10.098 


0.014 


0.00027 


11 


11.0223 


10.8 


0.0011 



In [IT] it was proved that the set of eigenvalues of the problem tends to 
Z as e ^ 0. Figure 1 (where the eigenvalues have been calculated by SPPS) 
illustrates this assertion. 

6 Conclusions 

In the present work a new representation for solutions of the Sturm-Liouville 
equation is presented, which we call a spectral parameter power series (SPPS). 
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0.5 0.4 0.3 0.2 0.1 0.05 0.03 

Figure 1: Convergence of A„ to n as e — > 0. 



It gives a new natural and powerful method for solving initial value, boundary 
value, and spectral problems. We have shown that it is applicable not only to 
regular problems but also to singular problems as well as to Sturm-Liouville 
problems with spectral parameter dependent boundary conditions. The nu- 
merical realization of the method is simple, and its ideas can be explained to 
undergraduate students in engineering, mathematics, and physics. One of the 
important advantages of the SPPS method consists in the automatic calculation 
of eigenfunctions together with the eigenvalues of a spectral problem. The SPPS 
method is competitive with the best software for numerical solution of ordinary 
differential equations and spectral problems currently available. The authors' 
goal has been to show the broad applicability and competitiveness of the method 
rather than to achieve spectacular numerical results; we believe that much bet- 
ter performance of the SPPS method can be achieved with the implementation 
of additional computational techniques, some of which have been mentioned in 
this paper. 

Acknowledgements: The authors express their gratitude to CONACYT 
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